####################################################################
# File-Name: ItPSR simulations for graph (2)                       #
# Date: November 2015                                              #
# Author: Ruth Dassonneville                                       #                      
####################################################################

## load packages
install.packages("foreign")
install.packages("plm")
install.packages("lme4")

library(foreign)
library(plm)
library(lme4)


## Logit model explaining voting for incumbent
## read data
data <- read.dta("/Users/ruthdassonneville/Desktop/ItPSR_Punishing Local Incumbents_replication.dta")

## run model
M1 <- glmer(IncumbencyCoalition2 ~ LaggedVote + Age+ Female + French +
              HighEdu+ Years_municipality+ Left_right + Terms + ENP + Unempl_change_10_11
            + HighEduUnempl 
            + (1+HighEdu|Municipality_code), family="binomial", data=data) 
summary(M1)

##get coefficients and varcov matrix
est<-fixef(M1)
V<-vcov(M1)
est
V

## set up sampling distribution
library(MASS)
N<-10000
S<-mvrnorm(10000, est,V)
head(S)

##choose sensible covariates
means<-apply(model.matrix(M1), 2, mean)
means

## predictions Unempl [11] -2 to 10 [adjust several times for each point prediction; first range for HighEdu=0, subsequently range for HighEdu=1]
## HighEdu [6] 0 or 1 [first series with HighEdu=0, subsequently adjust to HighEdu=1]
## Interaction [12] [adjust several times for each point predictions, range for HighEdu=0 and range for HighEdu=1]

X2<-means
X2
X2[11]<-10
X2[6]<-1
X2[12]<-10

## get linear predictor and form expected value
mu<-S%*%X2
mean(mu)
hist(mu)

lower<-apply(mu,2,quantile, 0.20)
lower
upper<-apply(mu,2,quantile, 0.80)
upper

#in response function to get predicted probability
p<- 1/(1+exp(-mu))
mean(p)

plower<-1/(1+exp(-lower))
pupper<-1/(1+exp(-upper))
plower
pupper

